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ABSTRACT 

We develop and test a method to compute mass and auto-correlation functions of rich clusters of 
galaxies from linear density fluctuations, based on the formalism of Gaussian peaks (Bardeen et al. 
1986). The essential, new ingredient in the current approach is a simultaneous and unique fixture of the 
size of the smoothing window for the density field, r/, and the critical height for collapse of a density 
peak, 5 C , for a given cluster mass (enclosed within the sphere of a given radius rather than the virial 
radius, which is hard to measure observationally). Of these two parameters, r/ depends on both the mass 
of the cluster in question and f2, whereas 5 C is a function of only $1 and A. These two parameters have 
formerly been treated as adjustable and approximate parameters. Thus, for the first time, the Gaussian 
Peak Method (GPM) becomes unambiguous, and more importantly, accurate, as is shown here. 

We apply this method to constrain all variants of the Gaussian cold dark matter (CDM) cosmological 
model using the observed abundance of local rich clusters of galaxies and the microwave background 
temperature fluctuations observed by COBE. The combined constraint fixes the power spectrum of any 
model to ~ 10% accuracy in both the shape and overall amplitude. To set the context for analyzing CDM 
models, we choose six representative models of current interest, including an J7q = 1 tilted cold dark 
matter model, a mixed hot and cold dark matter model with 20% mass in neutrinos, two lower density 
open models with £lo = 0.25 and Slo = 0.40, and two lower density flat models with (f2n = 0.25, Ao = 0.75) 
and (rio = 0.40, Ao = 0.60). This suite of CDM models should bracket any CDM model that is currently 
viable. The parameters of all these models are also consistent with a set of other constraints, including 
the Hubble constant, the age of the universe and the light-element nucleosynthesis with f2& chosen to 
maximize the viability of each model with respect to the observed gas fraction in X-ray clusters. 

Subject headings: Cosmology: large-scale structure of Universe - cosmology: theory - galaxies: 
clustering - galaxies: formation - numerical method 



1. INTRODUCTION 

Clusters of galaxies are cosmologically important be- 
cause they contain vitally important information on scales 
from a few megaparsecs to several hundred megaparsecs, 
and provide fossil evidence for some of the basic cosmolog- 
ical parameters (Richstone, Loeb, & Turner 1992; Bahcall, 
Fan, & Cen 1997). The fact that they are among the most 
luminous objects in the universe renders them an effec- 
tive and economical tracer of the large-scale structure, not 
only of the local universe (Bahcall 1988) but also of the 
universe at moderate-to- high redshift. 

Perhaps more interesting is the fact that clusters of 
galaxies are intrinsically rare with typical separations of 
~ 50/i _1 Mpc at z ~ and seemingly rarer at high red- 
shift (Luppino & Gioia 1995; Carlberg et al. 1996; Post- 
man et al. 1996). Their rarity is traceable to the fact that 
they only form at the rare, high peaks in the initial den- 
sity field. Since the mass in a sphere of radius 10/i _1 Mpc 
roughly corresponds to the mass of a rich cluster like the 
Coma cluster, the abundance of clusters of galaxies (i.e., 
the cluster mass function, Bahcall & Cen 1992) provides 
a sensitive test of the amplitude of the density fluctua- 
tions on that scale and places one of the most stringent 
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The spatial distribution of clusters of galaxies provides 
complementary information for cosmological models. A 
widely used statistic for clusters of galaxies is the two- 
point auto-correlation function. Earlier pioneering work 
(Bahcall & Soneira 1983; Klypin & Kopylov 1983) has 
met with dramatic improvements in recent years thanks 
to larger and/or new cluster samples that have become 
available (Postman, Huchra, & Geller 1992; Nichol et al. 
1992; Dalton et al. 1992,1994; Romer et al. 1994; Croft et 
al. 1997 ). Comparing with cosmological models clearly 
show that the two-point correlation function of clusters of 
galaxies provides a strong test on cosmological models on 
scales from several tens to several hundred megaparsecs 
QBardccn, Bond, fc Efstathiou 1987| ; Bahcall & Cen 1992; 
Mann, Heavens, & Peacock 1993; Holtzman & Primack 
1993; Croft & Efstathiou 1994; Borgani et al. 1995). In 
addition, the recent studies of superclusters and super- 
voids by Einasto et al. (1997a, b,c) show a very intriguing 
property that the correlation function of rich clusters ap- 
pears to be oscillatory on large scales. If confirmed, this 
would challenge most models. 

Hence, the combination of cluster mass and correla- 
tion functions provides a critical constraint on cosmolog- 
ical models on scales > 10 /i _1 Mpc. While uncertainties 
remain in the current clustering analyses as well as the 
abundance of observed clusters due chiefly to still limited 
cluster sample sizes typically of order of a few hundred 
clusters, large surveys underway such as the Sloan Digital 
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Sky Survey (SDSS; Knapp, Lupton, & Strauss 1996) and 
2dF galaxy redshift survey (Collcss 1998) should provide 
much more accurate determinations of both. 

The groundwork for the gravitational instability picture 
of cluster formation was laid down more than two decades 
ago (Gunn & Gott 1972). In the context of Gaussian cos- 
mological models, Kaiser (1984), in a classic paper, put 
forth the "biased" structure formation mechanism, where 
clusters of galaxies were predicted (correctly) to form at 
high peaks of the density field to explain the enhanced 
correlation of Abell clusters over that of galaxies. This 
idea was subsequently extended to objects on other scales 
including galaxies and the properties of linear Gaussian 
density fields were worked out in exquisite detail (Peacock 
& Heavens 1985; Bardeen et al. 1986, BBKS hereafter). 

While alternatives exist (e.g., Zel'dovich 1980; Vilenkin 
1981,1985; Turok 1989; Barriola & Vilenkin 1989; Bennett 
& Rhie 1990), a Gaussian model is simple and attractive 
(largely because of it) in that all its properties can be fully 
specified by one single function, the power spectrum of its 
density fluctuations. Moreover, it is predicted that ran- 
dom quantum fluctuations generated in the early universe 
naturally produced Gaussian density fluctuations, whose 
scales were then stretched to the scales of cosmological 
interest by inflation (Guth & Pi 1982; Albrecht & Stein- 
hardt 1982; Linde 1982; Bardeen, Steinhardt & Turner 
1983). Furthermore, observations of large-scale structure 
and microwave background fluctuations appear to favor a 
Gaussian picture (Vogclcy et al. 1994; Baugh, Gaztanaga, 
& Efstathiou 1995; Kogut et al. 1996; Colley, Gott, & Park 
1996; Protogeros & Weinberg 1997; Colley 1997). So mo- 
tivated, the present study will focus on the family of Gaus- 
sian CDM models. The reader is referred to Cen (1997c) 
for a discussion of the cluster correlation function in non- 
Gaussian models. We will employ the formalism of BBKS 
of Gaussian density field to devise an analytic method that 
can be used to directly compute the mass and correlation 
functions of clusters of galaxies. The needed input are: Pk 
(the power spectrum), flo and Ao. The method developed 
is calibrated and its accuracy checked by a large set of 
N-body simulations. 

The motivation for having such an analytic method is 
not only of an economical consideration (fast speed and 
much larger parameter space coverage possible) but also 
a necessity, especially for studying very rich clusters. For 
example, for clusters of mean separation of 200/i~ 1 Mpc 
(about richness 3 and above; Bahcall & Cen 1993, BC 
henceforth), a simulation box of size 1170/i _1 Mpc on a side 
would contain 200 such clusters, a number which may be 
required for reasonably sound statistical calculations. As- 
suming that the mass of such a cluster is 1.0 x 10 15 /i _1 M Q 
(approximately the mass of a richness 3 cluster; BC) and 
one requires 500 particles to claim an adequate resolu- 
tion of the cluster, it demands a requisite particle mass of 
2.0 x 10 12 /i _1 M Q . This particle mass requirement dictates 
that one discretize the whole simulation box into 10 8 ' 3 f2o 
particles (fto is the density parameter of a model) . Mean- 
time, a minimum nominal spatial resolution of 0.5/i _1 Mpc 
is needed to properly compute just the cluster masss within 
Abell radius of 1.5/i _1 Mpc, which translates to a spatial 
dynamic range of 2340. The combined spatial and mass 
resolution requirements are formidable for either PM code, 
which requires a very large mesh (2340 3 ) thus needs more 



than 57GB of RAM to allow for such a large simulation 
and hence is very expensive, if possible, or an adaptive 
code such as P 3 M (Efstathiou et al. 1985) or TPM code 
(Xu 1995), where CPU cost will be prohibitively large even 
if RAM permits. 

The paper is organized as follows. Descriptions of GPM 
for computing cluster mass function are presented in §2.1. 
Descriptions of GPM for computing the cluster corre- 
lation function are presented in §2.2. A calibration of 
Press-Schechter method using the fitted GPM parameters 
and some comparisons between GPM and Press-Schechter 
method are presented in §2.3. We discuss the various fac- 
tors that affect the cluster mass function in §3. Detailed 
constraints by the local rich clusters and the COBE obser- 
vations (Smoot et al. 1992) on all CDM models are pre- 
sented in §4. A simple as — relation (with errorbars) 
for CDM models is presented in §5, derived from fitting to 
the observed local cluster abundance alone. Conclusions 
are given in §6. 

2. GAUSSIAN PEAK METHOD FOR CLUSTERS OF 
GALAXIES 

2.1. Gaussian Peak Method for Cluster Mass Function 

It is convenient to define some frequently used symbols 
first. Hubble constant is H = 100/ikm/s/Mpc. f^o and Ao 
are the density parameters due to non-relativistic mater 
and cosmological constant, respectively, at redshift z = 0. 
fl z and A 2 are the same parameters at redshift z. r a is the 
comoving radius of a sphere in units of /i _1 Mpc, which in 
most times represents the Abell radius with value 1.5. r v is 
the virial radius in comoving /i _1 Mpc. r/ is the radius of a 
smoothing window in comoving /i~ 1 Mpc. Ma is the mass 
within a sphere of radius ta in units of h~ x M . M v is the 
mass within a sphere of radius r v (virial mass of a halo) in 
units of h~ 1 M . For formulae related to Gaussian density 
field we will follow the notation of BBKS throughout this 
paper. 

The cluster mass function may be derived by relating 
the initial density peaks to the final collapsed clusters, 
provided that peaks do not merge. Two pieces of observa- 
tions suggest that merger of initial density peaks of cluster 
size be infrequent. First, the typical separation between 
clusters of galaxies is ~ 100 h~ 1 Mpc, while the typical 
size of a cluster is ~ 1 /i Mpc. Second, empirical ev- 
idence of matter fluctuations, as indicated by observed 
galaxy number fluctuations (Davis & Peebles 1983; Strauss 
& Willick 1995), suggests that the current nonlinear scale 
is ~ 8 /i~ 1 Mpc, which is just about the size of fluctuations 
that collapse to form clusters of galaxies; i.e., the majority 
of clusters of galaxies form at low redshift. However, a 
more quantitative argument, that merger rate should be 
small, can be made as follows. Suppose that a cluster is 
moving at velocity vo at z = 0, then we can compute the 
total comoving distance that the cluster has travelled in 
its entire lifetime as 

d cm = / Vof(t)(l + z)dt, (1) 

Jo 

where to is the current age of the universe, f(t) is a func- 
tion to describe the evolution of the (proper) peculiar ve- 
locity of the cluster and the last term (1 + z) relates co- 
moving distance to proper distance. To illustrate the point 
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we will first use fio = 1, which gives the following simple 
relation: 

d cm = v H Q . (2) 

To arrive to the above relation we have used the follow- 
ing simple relations: t = t (l + z)~ 3/2 , t = 2H ~ 1 /3, 
f{z) = (1 + z)~ x / 2 (linear growth rate of proper peculiar 
velocities; Peebles 1980). For any reasonable model clus- 
ters of galaxies do not move at a speed (peculiar velocity) 
much higher than <~ lOOOkm/s at present (Cen, Bahcall & 
Gramann 1994); it can be obtained approximately in lin- 
ear theory by integrating a power spectrum, smoothed by 
an appropriate window, to yield the total kinetic energy 
(e.g., Suto, Cen, & Ostriker 1992). Note that some galaxy 
in a virialized cluster may move at a higher speed, but 
we are not considering such objects. So, a cluster mov- 
ing at 1000 km/s today moves a total comoving distance 
of 10ft- _1 Mpc in an J7o = 1 universe. The same cluster 
will move a longer distance in a low flo universe, but not 
by a large factor. An upper bound on d cm in such cases 
may be obtained by setting £1 = 0, in which case we have 
f(z) = (1 + z), t Q = Hq 1 and t = t (l + z) _1 . The upper 
bound is 

dcm.ub = VoH (l + Z max ), (3) 

where z max is the maximum redshift to which O = ap- 
plies. Let us make a simple, approximate estimate for a 
realistic lower bound by taking £1 = 0.2, as follows. For 
an fio = 0.2 model, the redshift at which il = 0.5 is 3.0, 
which we denote as z max . We treat the redshift range 
z > z max as an fl = 1 model and treat z < z max as an 
£1 = model, which gives a simple result: 

d cm = v H [(1 + -niu.r 

) + !]■ (4) 

For z max = 3.0, d cm = BvqHq. Since velocity decayed 
from z max to z — 0, a more reasonable upper bound on 
v is 1000/(1 + z max )km/s. This gives 12.5ft, _1 Mpc for 
vq = lOOOkm/s and z max = 3. 

Since one needs to collapse a sphere of 9.5f2 1//3 /i _1 Mpc 
in a uniform density field to form a massive cluster of mass 
1 x 10 15 /i _1 M Q , i.e., cluster density peaks have to have a 
separation of at least <~ 20£1 /i _1 Mpc and more likely 
~ 50/i _1 Mpc (mean separation of rich clusters today), it 
thus seems quite unlikely that a significant fraction of any 
massive cluster peaks have merged by z = 0. This con- 
clusion is, however, not in conflict with observations that 
seem to show signs of recent and/or ongoing merger ac- 
tivity. In general, merger is an ever-going processes (at 
least in the past) in any plausible (i.e., a plausible range 
in f2g) hierarchical structure formation model. But, these 
mergers or substructures seen in some clusters are sub rich 
cluster scale mergers, i.e., sub-peaks within a large cluster 
scale peak are in the process of merging, a result which 
is in fact expected if clusters have been forming in the 
recent past in a hierarchical fashion. To our knowledge, 
there is no major merger event of two massive clusters ob- 
served. For example, in 55 Abell clusters catalogued by 
Dressier (1980), there is no case of two massive clusters in 
the process of imminent merging, although there does seem 
to have significant substructures in a significant fraction 
of clusters in various optical studies (e.g., Geller & Beers 
1982; Dressier & Shectman 1988; West, Oemler, & Dekel 
1988). X-ray observations also show a large fraction of 



clusters with substructures (see Forman & Jones 1994 for 
a review). That being said, one needs to be extra cautious 
in interpreting such sub cluster scale merger/substructure 
events, due to unavoidable projection contaminations (see 
Cen 1997 for a thorough discussion of projection effects). 

Having shown that merger should be infrequent, the key 
link then is to relate a density peak of height 

v = F/a (5) 

to the final mass of a cluster defined within a fixed ra- 
dius, say, the Abell radius ta- Here, F is a density field 
smoothed by a window of size 77 and Co is the rms fluctua- 
tion of F. Gaussian smoothing window (in Fourier space) 

W(kr f ) = cxp(-r 2 k 2 /2), (6) 

will be used throughout this paper, because it guaran- 
tees convergence of any spectral moment integral with 
any plausible power spectrum. Top-hat smoothing does 
not have this feature. For the sake of definiteness and 
convenience in comparing with observations, we define a 
cluster mass, Ma, as the mass in a sphere of comoving 
Abell radius, ta — 1-5 /i _1 Mpc, in most cases. Cluster 
mass defined otherwise will be noted in due course. But 
the formalism developed here should be applicable for any 
plausible radius. 

For a spherical perturbation, the mean density within 
the virial radius (at redshift z) in units of the global mean 
density (at redshift z) can be parameterized by 

p v {£l Zl A z ) = 178fL7°- 57 C(f! z ,A z ). (7) 

For Q. z = 1 and A z = 0, it is well known that p v = 178, 
a result first derived by Gunn & Gott (1972). In equation 
(7) C(fl z , A z ) (a function of both Q z and A z ) has a value 
close to unity. It has been shown that C = 1 is a good 
approximation for both A z = model (e.g., Lacey & Cole 
1993) and A z + Q 2 = 1 model (e.g., ECF) for the range 
of fi z of interest (0.1 < 2 < 1.0). The mass within the 
comoving virial radius, r v , at redshift z, is therefore 

47T 

M v (zf~-rlp v (n z ,A z )p c n 

=2.058 x W 14 n z - 57 il rlC(n z ,A z ) (8) 

where p c is the critical density at z — 0. Since the spher- 
ical perturbation smoothed by a Gaussian window of co- 
moving radius 77 is, by assumption, just virialized at the 
redshift in question, we have another expression for M v as 
a function of 77: 

M v (z) = {2ir) 3/2 r 3 fPc n a = AMI x 10 12 O r) (9) 

Next, we need to relate M v to Ma, the latter of which is 
observationally more obtainable. We assume the following 
simple relationship: 

M A = M v ( r -^f- a . (10) 

This relation holds exactly, if the density profile of the 
cluster has the power-law form: 

p(r)ocr- a . (11) 
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But, in general, the density profile of a cluster does not 
have a power-law form, so a(Q z , A z , Ma, P^) should only 
be considered as a fitting parameter, which should, in prin- 
ciple, be dependent on both the cluster mass and under- 
lying cosmology. However, motivated by the insight of 
Navarro, Frenk, & White (1996) that there seems to be a 
universal function (as a function of scaled radius in units 
of the virial radius each individual halo) for density pro- 
files of dark matter halos, independent of cosmology and 
halo mass, it is hoped that a will only be a weak func- 
tion of both the underlying cosmology and cluster mass. 
As a matter of fact, as we will show below, the best fit 
to N-body results requires that a be a constant equal to 
2.3, in harmony with the work of Navarro ct al. (1996). 
Combining equations 8,9,10 yields 



r f ^17r^- 3)/a 



M A 



2.058 x 10 14 



l/a 



ft0.19(a-3)/a n -l/a C(f ^ ^(a-SJ/Sa (12) 

This equation allows us to determine the smoothing radius 
rf for a cluster of mass Ma (within radius r A ), given a. 
We see that r/ will be completely deterministic, if a and 
C can be specified a priori. As we will subsequently show, 
we will choose to fix C = 1 and let its dependence on f2 z 
and A z be absorbed by another fitting parameter S c (see 
below). Therefore, the final equation for 77 that we will 
use is 



r^l7r^- 3)/a 

n 0.19(a-3)/a fi - 



M A 



2.058 x 10 14 

l/a 



l/a 



(13) 



In this equation a is the only adjustable parameter. But as 
will be shown later, a turns out to be a constant. There- 
fore, 77 is no longer an adjustable parameter, rather it is 
a unique function of only Ma and rA ■ Another useful ex- 
pression is to relate M v to Ma (obtained by combining 
equations 9 and 13) in terms of rA- 



or 



M£&58 x 10 14 r^ (a - 3)/a 

^0.57(a-3)/Q^(a-3)/Q 



M A 2£58 x 10 14 r^ a 



M A 



2.058 x 10 14 



3/a 



(14) 



M v 



2.058 x 10 14 



a/3 



Q0.19(3-a)Q(3-a)/3 



(15) 



We note that, in the special case where rA — r v , we have 
(from equation 9) 



w > A 1/3 0-1/3 

Sin 



s 4.347xl0 12 ' 
which is independent of a, and the virial radius is 

1/3 



(16) 



M A 



2.058 x 10 14 



n°-"V 8 . (17) 



The circular velocity of the (just virialized) halo at redshift 
z can be expressed as 



GM v (l + z) 



1 1/2 



=£4 x 10 2 



M A 



1/3 



1.0 x 10 14 , 

O-0.095ol/6 (1 + z) l/2 km/ ^ 



(18) 



Note that the dependences of v c on f2 z and ilo are very 
weak. But v c depends rather strongly on z and somewhat 
strongly on M v . The 1-d velocity dispersion a\\ is just 
equal to v c /\/2. Another useful relation is between 77 and 
v c (or o-||): 



77 = 3.845 x I0~ 3 v c tt 



0.09.-, 



-1/2 



+ 



-1/2 



(19) 



where v c is in km/s. 

Once we have uniquely determined the smoothing win- 
dow for a given cluster mass and cosmology (equation 13), 
there is only one last parameter to be specified before the 
abundance of peaks forming clusters of a given mass is 
uniquely fixed: we need to specify the required peak height 
such that density peaks with such a height just collapse 
and virialize at the redshift in question. Here once again, 
we are guided by the spherical perturbation model, and use 
S c as a parameter to quantify the required peak height. 5 C 
is the linear overdensity of a peak at the concerned red- 
shift at which a peak just collapses and virializes. S c is 
1.68 in the spherical top- hat collapse in an Q z = 1 model 
(Gunn & Gott 1972). Since realistic density perturbations 
are likely to be non-spherical and non-top-hat and there 
are models other than those with Q z = 1 also being con- 
sidered, S c can only be treated as a fitting parameter, to 
be determined by comparing to N-body simulations. In 
other words, the spherical collapse model is not applied 
explicitly, rather it is used as a guide for an initial guess of 
S c . The actual collapse, or some portion of the collapse in 
time or in space, of a density peak may not be spherical. 
For example, a triaxial proto-cluster will require a lower 
value of S c than a spherical proto-cluster (More, Heavens, 
& Peacock 1986). The validity of our method does not 
critically depend on the spherical collapse model, rather it 
is made valid by comparing to N-body simulations treat- 
ing S c as a fitting parameter. Although S c could depend 
on Pk and cluster virial masss M v as well as fi z and A z , 
we will restrict our fitting procedure as if S c only depends 
on f2 z and A z . In any case, the final fit to N-body results 
proves that this assumption is good. 

We now set down the basic formulae from BBKS for 
computing the number density of peaks of appropriate 
sizes and heights for a Gaussian density field. The dif- 
ferential peak density is 



N pk (u) 



1 



(27T) 2 i?3 



- 2 / 2 G( 7 ,7^), 



(20) 



where G(7, w) is, for the convenience of calculation, an an- 
alytic formula approximating the exact three-dimensional 
integral (see equation A12 of BBKS). G(-f,w) is accurate 
to better than 1% over the range 0.3 < 7 < 0.7 and 
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— 1 < w < oo with the accuracy being better than 0.1% 
for w > 1, according to BBKS: 



_ w 3 - 3 7 2 w + [5(7)w 2 + Ci(-y)] cxp[-^( 7 )w 2 ] 
(7 ' Wj ~ 1 + C 2 ( 7 )exp[-C 3 (7H 

(21) 

(equation 4.4 of BBKS), where to = ju. For all the models 
which we have computed or are of current interest we have 
7 > 0.7 and w > 1.1 in most cases with the lowest values 
being 7 = 0.58 and w = 0.84. Therefore, it is accurate 
to use G(j, to) for models of current interest. The various 
symbols arc defined as follows. 



A = 



D 



5/2 



9-5 7 2 



432 



(107r) 1 /2(9_ 5 7 2)5/2 

Ci = 1.84+ 1.13(1 - 7 2 ) 5 - 72 
C* 2 =8.91 + 1.27exp(6.5l7 2 ) 
C 3 = 2.58cxp(1.057 2 ) 



(22) 



(equation 4.5 of BBKS). The parameters 7 and i?* are 
related to the moments of the power spectrum: 



7 o-i 



0-2 



(23) 



(equation 4.6a of BBKS), and Uj are spectral moments: 
*i ee / ^P kk ^W 2 (kr f ) (24) 

(equation 4.6c of BBKS). We are now ready to derive the 
cluster mass function by counting peaks of appropriate 
sizes and heights. The procedure can be described in four 
steps: 

1. For a cluster mass Ma, we find rj using equation 13. 

2. We smooth the power spectrum Pk by the square of a 
Gaussian window of radius 77 (equation 6) and then com- 
pute Co, 7 and i?* (equations 23,24). 

3. Requiring that vao = 5 C yields the threshold peak 
height v t = S c /a . 

4. Integrating equation 20 from v t to 00 gives the cumu- 
lative cluster mass function: 



/>oo 

n(> M A ) = / N pk (v)dv. 

J Vt 



(25) 



To summarize, we have three adjustable parameters: 
a{Q z ,A z ,M A ,Pk), 5 c (n z ,A z ) and C{Ct z ,A z ). In general, 
a(Q z , A z , Ma, Pk) should be a function of four variables, 
Q z , A z , Ma and Pk , whereas <5 C and C should depend only 
on Q z and A z . We choose to fix C(Cl z ,A z ) to be unity, 
independent of cosmological parameters, and treat only a 
and S c as two adjustable parameters. Since both parame- 
ters (C and 5 C ) depend on the same cosmological param- 
eters, this treatment is justified and simplifies the fitting 
procedure; the dependence of C(Q Z ,A Z ) on cosmological 
parameters is assumed to be absorbed by 5 C (Q Z ,A Z ). 



We now turn to N-body simulations to calibrate GPM to 
compute the cluster mass function, i.e., to determine the 
two fitting parameters in GPM, a and 5 C . We have thirty 
two cosmological models at our disposal to calibrate GPM 
and test its accuracy. The models are listed in Table 1. 
The second and third columns are the density parame- 
ter and cosmological constant of the model, respectively. 
The fourth column is the linear rms density fluctuation 
in an 8 /i _1 Mpc top-hat sphere at z = 0. The fifth col- 
umn, EP (Excess Power), is a parameter to describe the 
shape of the power spectrum on scales ~ 8 — 300 /i _1 Mpc, 
introduced by Wright et al. (1992). The definition is 
EP ee 3.4<725/(78, where 025 is the linear rms density fluc- 
tuation in a 25 /i _1 Mpc top-hat sphere at z = 0. The last 
column indicates the type of power spectrum used (details 
are given in the table footnotes). 

Some models are physically self-consistent in the sense 
that their power spectrum transfer functions are computed 
for the given cosmological parameters, while others are 
not. The latter are included to increase the coverage of the 
parameter space for calibration purposes. Taken together 
these thirty two models span the ranges of fl z , A z , cr 8 and 
EP of current interest: fl z = 0.2 -> 1.0, A z = 0.0 -> 0.8, 
<7 8 = 0.35 -> 1.5, and EP = 0.735 -> 1.923. As a note, 
Wright et al. (1992) find that the range of EP that fits 
the COBE data is 1.30 ± 0.15 (1<t), which is consistent 
with analysis of the galaxy power spectrum by Peacock 
& Dodds (1994) and Feldman, Kaiser & Peacock (1994), 
or observations of large-scale galaxy clustering by Maddox 
et al. (1990). The reason that we use EP rather than T 
(= Qoh) is that T can only be used for CDM type models 
with n = 1. 

Each model is run using a particle-mesh code with a 
box size 400 /i~ 1 Mpc. A large simulation box is needed 
in order to produce a significant number of the rich but 
rare clusters. The simulation box contains 720 3 cells and 
240 3 = 10 7 1 dark matter particles, with a particle mass of 
1.3 x 10 12 f2o h~ l Mq. The mass resolution seems adequate 
for our purpose: a cluster of mass 6 x 10 14 /i _1 M contains 
462£!q 1 particles. In each simulation, clusters are selected 
as the maxima of the mass distribution within spheres of 
comoving radius of ta = 1.5h _1 Mpc. The mass of each 
cluster is determined in a sphere within the Abell radius 
ta- The results are not sensitive to the cluster-finding 
algorithm as long as the mass is defined within a chosen 
radius. In other words, different group-finding algorithms 
such as friends-of-friends or DENMAX are all able to lo- 
cate the density peaks properly. 

Before we present our numerical results, it is important 
to check resolution effect. We make the following reso- 
lution calibration test. We run two simulations for the 
SCDM model (model 1 in Table 1) with an identical ini- 
tial realization in a box of size L = 128/i _1 Mpc. One of the 
two simulations uses the same PM code as for the thirty 
two models listed in Table 1 and has the same numer- 
ical resolution (1.39 /i~ 1 Mpc); the other simulation has 
a much higher resolution (0.0625 /i _1 Mpc) based on the 
P 3 M scheme utilizing a special computer chip (GRAPE) to 
solve the PP part of the force computation (Brieu, Sum- 
mers & Ostriker 1995). Since the resolution element of 
the P 3 M simulation is much smaller than the Abell ra- 
dius, it can be considered as having an infinite resolution 
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Table 1 

List of parameters for 32 models 



ivioaci 


o 


A 


°8 


hjF = ,5. 4(725 / 08 


Comment 


1 


1.000 


0.000 


1.050 


0.956 


SCDM a 


2 


1.000 


0.000 


0.700 


0.956 


SCDM a 


3 


1.000 


0.000 


0.525 


0.956 


SCDM a 


4 


1.000 


0.000 


0.350 


0.956 


SCDM a 


5 


1.000 


0.000 


1.050 


1.088 


Pk = k- 1 


6 


1.000 


0.000 


0.808 


1.088 


Pk = k- 1 


7 


1.000 


0.000 


0.700 


1.088 


Pk = k- 1 


8 


1.000 


0.000 


0.525 


1.088 


Pk = k- 1 


9 


1.000 


0.000 


0.350 


1.088 


Pk = k' 1 


10 


1.000 


0.000 


1.050 


1.923 


Pk = k' 2 


11 


1.000 


0.000 


0.808 


1.923 


Pk = k- 2 


12 


1.000 


0.000 


0.700 


1.923 


P k = k- 2 


13 


1.000 


0.000 


0.525 


1.923 


Pk = k- 2 


14 


0.350 


0.000 


0.800 


1.196 


OCDMl b 


15 


0.410 


0.000 


0.689 


1.196 


OCDM1& 


16 


0.446 


0.000 


0.632 


1.196 


OCDMl b 


17 


0.520 


0.000 


0.524 


1.196 


OCDMl b 


18 


0.600 


0.000 


1.000 


1.060 


OCDM2 b 


19 


0.661 


0.000 


0.818 


1.060 


OCDM2 b 


20 


0.692 


0.000 


0.730 


1.060 


OCDM2 c 


21 


0.750 


0.000 


0.575 


1.060 


OCDM2 c 


22 


0.818 


0.000 


0.404 


1.060 


OCDM2 c 


23 


0.400 


0.600 


0.790 


1.225 


LCDMl d 


24 


0.692 


0.308 


0.591 


1.225 


LCDMl d 


25 


0.842 


0.158 


0.460 


1.225 


LCDMl d 


26 


0.200 


0.800 


1.500 


1.225 


LCDM2° 


27 


0.355 


0.645 


1.321 


1.225 


LCDM2° 


28 


0.458 


0.542 


1.211 


1.225 


LCDM2 6 


29 


0.573 


0.427 


1.088 


1.225 


LCDM2° 


30 


0.667 


0.333 


0.982 


1.225 


LCDM2° 


31 


0.796 


0.204 


0.813 


1.225 


LCDM2 6 


32 


0.871 


0.129 


0.690 


1.225 


LCDM2 e 



a the standard CDM model with Hubble constant H Q = 
50km/s/Mpc, Oo = 1-0 and n = 1.0, where n is the power in- 
dex at very large scale. BBKS power spectrum transfer function 
(equation G3) is used 

b an open CDM model with Hubble constant H a = 70km/s/Mpc, 
O = 0.35 and n — 1.0; BBKS power spectrum transfer function 
(equation G3) is used. 

c an open CDM model with Hubble constant H a = 60km/s/Mpc, 
O = 0.60 and n = 1.0; BBKS power spectrum transfer function 
(equation G3) is used. 

d a CDM model with a cosmological constant with Hubble con- 
stant H = 65km/s/Mpc, fi = 0.40, A = 0.60 and n = 0.95; 
the power spectrum transfer function is computed as in Cen et al. 
1993. 

e a CDM model with a cosmological constant with Hubble con- 
stant H = lOOkm/s/Mpc, fi = 0.20, A = 0.80 and n = 0.95; 
the same power spectrum as LCDM1 is used. 
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Fig. 1. — Mhr/Mlr is plotted against M^r for the three models. The solid line in each plot is the best linear fit: Mur/M^r = 
a + blog 10 M LR , where M LR is in fc- 1 M . We find (a, 6) to be (-0.469,0.103), (-0.0528,0.0662), (-0.952,0.137) for models (1,14,29), 
respectively. 
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Fig. 2. — Mass functions of various models. Each curve is labeled by its model number from Table 1. The simulation results are shown as 
symbols with horizontal errorbars being the uncertainties in the mass determination (15%) and the vertical errorbars being the statistical la 
errorbars for the number of clusters. The solid curves in Figure 2 are the results from GPM. Note that simulation box size limits the density 
of clusters to > 1.56 X 10 — 8 /i 3 Mpc — 3 , when there is only one cluster in the whole box. At /(> Ma) = 10- 6 fe 3 Mpc- 3 , there are 65 clusters 
in the simulation box. 



(i.e., representing the truth) for the purpose of calibrating 
our low resolution results. Because the two simulations 
have identical initial conditions, we are able to identify 
every rich cluster in one simulation with its counterpart in 
the other. Having made such a one-to-one correspondence 
we can compute the ratio Mhr/Mlr (HR stands for high 
resolution and LR for low resolution) as a function of clus- 
ter mass Mifl. This allows us to make corrections to Mlh 
in the lower resolution simulation. The above resolution 
calibration procedure is repeated for an open CDM model 
(model 14 in Table 1) and a CDM model with a cosmo- 
logical constant (model 23 in Table 1). 

Figures (la,b,c) show the results for the three mod- 
els, where Mhr/Mlr is plotted against Mlr- The solid 
line in each plot is the best linear fit: Mhr/Mlr = 
a+61og 10 Mlr, where Mlr is in h~ x M Q . We find (a, b) to 
be (-0.469,0.103), (-0.0528,0.0662), (-0.952,0.137) for 
models (1,14,23), respectively. Let us call the three fitting 
functions (M hr /M L r) as R^M), R 035 (M), R 04 {M), for 
the three models run: Cl z = 1 (model 1 in Table 1), 2 = 
0.35 and A z = (model 14), and Q 2 = 0.40 and A z = 0.60 
(model 23). Then, for the mass of each cluster, Mlr, we 
correct it by multiplying it by i?i + (i? 35— -Ri)(l— ^ z )/0.65 
in an open model, or by R\ + (i?04 — i?i)(l — O z )/0.60 in a 
A model, where il z is the density parameter of the model 
under consideration. From Figure 1 we see that the typi- 
cal correction is about 5-10% in the upward direction with 
a dispersion of ~ 5%; i.e., lower resolution simulations 
slightly underestimate the mass within the Abell radius, 
as expected. Calibrating the lower resolution simulation 
results we assign an errorbar of 15% for each cluster mass. 



We are now ready to find the best fitting parameters 
(a, S c ) by comparing results from GPM to the direct N- 
body results. Before starting the fitting procedure, we 
have some rough idea about what the values of a and 5 C 
may be. We pick a = 2.5 and 6 C = 1.5 as an initial guess. 
In the end, the best values are found to be 

a = 2.3, (26) 

a constant independent of the cluster mass and cosmology, 
and 

x _/l-40- 0.01(1.0 -Sl x ) for A z =0 , , 

c \1.40 + 0.10(1.0 -fl z ) for Q z + A z = l. [ <) 

The best overall fit for all the models is judged by the au- 
thor by direct visual examination. We find it very difficult 
to design an automated fitting procedure to be gauged by 
some objective parameters, because of the enormous range 
of the number densities of clusters and hard-to-define er- 
rorbars hence weighting schemes for the densities. But as 
we will see, the final fits are probably as good as one would 
have hoped, which suggests that our somewhat subjective 
fitting procedure works very well. In any case, the final fit 
values span very narrow ranges (in fact, a turns out to a 
constant, and S c varies from 1.40 to 1.39 from fl z = 1 to 
in A z = models, from 1.40 to 1.36 from Cl z = 1 to in 
A z + A z = 1 models), which indicates that the fitting pro- 
cedure is robust and stable. We note that the fitted value 
of a is in excellent agreement with observations (Carlberg 
et al. 1996; Fischer et al. 1997). It seems useful to es- 
timate the uncertainties in the fitted values of a and S c . 
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However, the sensitivity of a fit to the two parameters de- 
pends on the fitted mass function itself: a low amplitude 
mass function depends more sensitively on the two param- 
eters than a high amplitude mass function. This is so, of 
course, because the abundance of rarer objects depends 
more sensitively on the parameters. Roughly speaking, 
a serves more to fix the shape of the mass function in a 
somewhat less sensitive way, while S c determines the over- 
all amplitude and more sensitively the amplitude on the 
high mass end of the mass function. Our estimates on the 
uncertainties are Aa = 0.1 and A5 C = 0.01. 

Figures (2a, b) show the simulation results as symbols 
for 32 models with horizontal errorbars being the uncer- 
tainties in the mass determination (15%) and the ver- 
tical errorbars being the statistical lcr errorbars. Note 
that simulation box size limits the density of clusters to 
> 1.56 x 10~ 8 /i 3 Mpc~ 3 , when there is only one cluster in 
the whole box. At /(> M A ) = 10" 6 ft, 3 Mpc" 3 , there are 65 
clusters. The solid curves in Figure 2 are the results from 
GPM. We see that the GPM results fit remarkably well the 
simulation results for all the thirty two models examined. 
We note that the actual errorbars should be larger than 
what are shown for the simulated results because of cosmic 
variances; i.e., the simulation boxsize, although quite large 
being 400 /i _1 Mpc, may still not be large enough to have 
the cosmic variance diminished, especially for models with 
significant power on several hundred megaparsecs scales. 
In any case, the GPM results fit the N-body results for all 
the models within 2a in the vertical axis, and within a fac- 
tor of 1.25 in the horizontal axis. Since the observed mass 
function (BC) has uncertainties in mass about a factor 2.0 
and in number density about a factor > 2.0, the GPM re- 
sults are practically precise, for the purpose of comparing 
model results computed using GPM with observations. 

At this point it seems appropriate to reiterate the virtue 
of the current method. The essential unique ingredients 
are the introduction of two adjustable parameters, a and 
8 C , the first of which turns out to be a constant and the 
second of which can be simply expressed as a function of 
Q z (the density parameter at the redshift in question). 
Note that the fitted parameter 6 c (Cl z ) has a slightly dif- 
ferent form for the case with A z = than for the case 
with fi z + A z = 1 . The fact that 8 C is only a rather weak 
function of fl z in both cases indicates that the method is 
robust. 

2.2. Gaussian Peak Method for Cluster Correlation 
Function 

Having found that the initial density peaks, appropri- 
ately defined, indeed correspond to the clusters formed at 
late times, as indicated by the goodness of the fits of the 
results from GPM to the N-body results in terms of clus- 
ter mass function presented in the preceding section, we 
have some confidence that we may be able to compute the 
cluster-cluster two-point correlation function using GPM. 
We will now proceed along this route. 

Even in the linear regime where dynamic contribution to 
the clustering can be ignored, the primary difficulty in cal- 
culating the cluster-cluster two-point correlation function 
using Gaussian peaks is the ambiguity of relating appro- 
priate peaks to the clusters of interest. This ambiguity is 
reflected in one's inability to fix the smoothing length rj 
and the threshold peak height v t . The GPM, described 



in §2, eliminates this ambiguity by simultaneously fixing 
both rj and v t . This is achieved by demanding that the 
appropriate peaks yield the correct cluster mass function, 
when compared to direct N-body simulations. 

Following BBKS, we use the following approximate for- 
mula, which is applicable when the correlation function is 
smaller than unity [however, BBKS state that it may well 
be a reasonable approximation even when the statistical 
correlation function (first term at the right hand side of 
equation 28, see below) is not really small], to compute 
the final cluster-cluster correlation function including lin- 
ear dynamical contributions: 

S, P k,pk ~ ^ < o _ o > + 1^ £p, P , (28) 

(equation 6.63 of BBKS), where £ Pt/) is the two-point den- 
sity auto-correlation function, and < v > is defined by 

/■OO 

v= / \v- 70/(1 -^Npkdv (29) 

J v t 

(equation 6.45 of BBKS), where v t is determined in the 
procedure described in §2; 7 is defined in equation 23; 
N p k is defined in equation 20. Note that we have chosen 
a step function for threshold function t{v/v t ) in equation 
6.45 of BBKS to arrive at equation 29, i.e., only peaks 
above v t are assumed to be able to collapse and no peaks 
below v t are allowed to collapse. A smoother threshold 
function (e.g., Kaiser & Davis 1985; BBKS) may be used, 
but the primary effect will be to slightly change the fitting 
parameter S c . We therefore will use the sharp step func- 
tion as the threshold function without loss of generality. 9 
is defined by 

a _ 3(1 - 7 2 ) + (1-216 - 0.9 7 4 ) exp[-7/2( 7 ^/2) 2 ] 
[3(1 - 7 2 )+ 0.45 +( 7 zV2) 2 ] 1 /2 +7 j,/ 2 

(equation 6.14 of BBKS). We note that equation 28 is 
valid in the linear regime when £ P:P is much less than 
unity and it is not yet clear whether the approximation, 
coupled with our definitive peak identification method, 
also works in the mild nonlinear regime. Our goal is to 
find an approximation based on equation 28 which will 
give sufficiently accurate results for £ p k, P k in the regime 
whose values are of order unity and below. For this reason 
we choose to modify equation 28 in the following man- 
ner: £ pk , p k ~ (^T + l) £ P , P D(£p,p), where D(£ PtP ) is 

a fitting parameter that depends only on £ PtP , which is 
the mean two-point matter correlation function, defined 
as £ p , p (x) = J* Z p ,p(y)y 2 dy. The reason for using £ PiP 
instead of £ PjP is that £ PiP is a better indicator of nonlin- 
earity than £ PjP . The form of D should be constrained at 
the linear end: D(0) = 1. As we will show below, it turns 
out that equation 28 fits results very well; i.e., fitting to 
numerical results indicates that D = 1 is a good approx- 
imation for the interested range in ^ p k, P k- To be clear, 
we use equation 28 for all the subsequent calculations of 
cluster correlation functions. 

We now compare the two-point peak-peak correla- 
tion function, calculated using GPM described above, to 
the cluster-cluster two-point correlation function obtained 
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Fig. 3. — The cluster correlation functions for clusters with mean separation of 55/i _1 Mpc, for 28 models as indicated in the panels. The 
errorbars are la statistical. Three curves are shown in each panel for each model; the middle curve is obtained using equation 28 and the top 
and bottom curves are obtained by adding ±2/i -1 Mpc to each point of the middle curve in the x-axis. 
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from N-body simulations. We compute the two-point cor- 
relation function from N-body simulations using the fol- 
lowing estimator: 



Ccc(r) 



N CR (r) 
N RR (r) 



1. 



(31) 



where Nc R (r) and Nnn(r) are the number of pairs be- 
tween clusters and random spatial points and the num- 
ber of pairs among random spatial points, respectively, 
at separation r — > r + Ar. The number of random spa- 
tial points for each realization within the simulation box 
(400/i _1 Mpc) is chosen to be the same as the number of 
clusters in question. A total of 100 random realizations are 
made. The final £ cc is averaged over the 100 estimations 
and errorbars are estimated by separately computing the 
correlations for each of the eight octants of each simulation 
box. 

Figures (3a,b,c,d,e,f,g) show the results of 28 models for 
clusters with mean separation of 55/i _1 Mpc. Three of the 
remaining four models do not have enough clusters whose 
masses exceed 1.4 x 10 14 /i _1 M Q (the low mass cutoff for 
clusters found in the simulations). Model 32 is not shown 
only to save space, although the goodness of the fit of GPM 
result to the N-body result in Model 32 is comparable to 
that of Model 15 (top left panel of Figure 3d). Figures 
(4a,b,c,d,e) show the results of 20 models for clusters with 
mean separation of 34ft,~ 1 Mpc. The models which are not 
shown, again, do not have enough clusters with masses 
greater than 1.4 x lO 14 /*" 1 M Q , except for Model 32. The 
errorbars are la statistical. Three curves are shown in each 
panel for each model; the middle curve is what is obtained 
using equation 28 and the top and bottom curves are ob- 
tained by adding ±2/i~ 1 Mpc to each point of the middle 
curve in the x-axis. We see that GPM works well in the 
range of scales where £ cc = 0.1 — 2.0 for all the models. 
Some models fit for even larger ranges. In particular, the 
correlation length seems quite accurately computable by 
GPM; r (GPM)±2h- 1 Mpc [where r (GPM) is the length 
scale where the correlation computed by GPM (equation 
28) is unity] appears in agreement with the correlation 
length computed from direct N-body simulations for all 
the models studied here. 

Several authors have used peaks in Gaussian density 
fields to compute the cluster-cluster two-point correlation 
function. Mann, Heavens, & Peacock (1993) use a method 
developed by Bond & Couchman (1988) which combines 
the Gaussian peak formalism with Zel'dovich approxima- 
tion. Their method is analytically tractable and fast. Un- 
fortunately, the accuracy of their method has not been 
carefully checked by N-body simulations and the two pa- 
rameters rf and S c (or r/ and v t ) are not fully determin- 
istic. It will be useful to apply the treatment of r/, a and 
S c here to their method. Holtzman & Primack (1993) use 
a similar but somewhat different formula (Bardeen, Bond, 
& Efstathiou 1987) than what is used here to compute 
the cluster correlation function in some variants of CDM 
models. Aside from the slightly different formula used to 
compute the statistical and dynamical correlation terms, 
the primary difference between ours and theirs lies in the 
treatment of r/ and 5 C . In their work rf is determined by 
the mass of a cluster in a rough way without taking into 
account the fact that the virial radius of a collapsed object 



is different from the radius within which observed mass is 
defined, which consequently also affects 5 C . It might be 
beneficial to check their analytic method against N-body 
simulations for a wide range of models, before firm con- 
clusions from detailed comparisons between models and 
observations can be drawn. Croft & Efstathiou (1994) 
have checked the same formula as used here to compute 
the correlation function of clusters using N-body simula- 
tions. The difference between our method and theirs is 
that r/ and 5 C in their approach are only roughly guessed 
based on a consideration of the approximate mass of clus- 
ters which determines 77 and the spherical collapse model 
which determines S c = 1.69. Their Gaussian peak method 
agrees well with their N-body results for the poorer clus- 
ters (d c = 15ft. _1 Mpc) in all three models (f2g = 1 model, 
£!o = 0.2 and A = 0.8 model, and f^o = 0.2 and A model) 
but the correlation functions from their Gaussian method 
are significantly lower than those from N-body simulations 
for richer clusters (d c — 35, 50/i _1 Mpc). The behavior of 
their results can be explained by the fact that they choose 
rf to be independent of Ma (or d c , the mean separation 
of clusters under consideration; a larger d c corresponds to 
richer, more massive clusters), whereas in our case r/ cor- 
relates with Ma- Consequently, the correlations of richer 
clusters from the Gaussian method of Croft & Efstathiou 
(1994) are underestimated. The primary drawback of all 
these previous studies is the inability to uniquely identify 
a set of peaks with a set of clusters of interest. Again, 
this traces to the ambiguity of choosing r/ and v t . Our 
method eliminates this ambiguity entirely and the results 
seem satisfactory. 

2.3. Calibrating Press-Schechter Formalism 

We have shown in §2.1 and §2.2 that GPM works quite 
well for computing both the cluster mass and correlation 
functions. It is tempting to try the same arguments on the 
widely used Press-Schechter (1974; PS) formalism. We re- 
call that the essential ingredient in GPM is the introduc- 
tion of two parameters: a and S c . Fitting results from 
GPM to N-body results fixes a to be a constant of 2.3 and 
6 C as a function of fl z (equation 27). The fitted value of 
a is consistent with both the theoretical work by Navarro, 
Frcnk & White (1996) and observations (Carlberg et al. 
1996; Fischer et al. 1997). So, we will adopt the same a 
derived from GPM to calibrate the PS formalism. But we 
will adjust 5 C by comparing PS results to N-body results. 

The basic PS ansatz results in the differential virialized 
halo mass function as: 

,,, U w /2 P S c da (M v ) 
n{M v )dM v -^j 3 — —27 — — 



exp(- 



7r M v o~q(M v ) dM v 
81 



2a 2 (M v ) 



)dM v , 



(32) 



where p is the mean density of the universe at the red- 
shift under consideration. Substituting M v by Ma using 
equation 14 with a — 2.3 gives 

n(M A )dM A -^.M5 x lO 15 ^ 913 ^ 173 ^ 304 

-1.304 8 C dao(MA) 



M 



exp( 



o-Z(Ma) dM A 



2al{M A ) 



)dM A - 



(33) 
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Fig. 5. — Mass functions of various models. The simulation results are shown as symbols with horizontal errorbars being the uncertainties 
in the mass determination (15%) and the vertical errorbars being the statistical la errorbars for the number of clusters. The solid curves are 
the results from Prcss-Schcchtcr method. 
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To calculate (Tq(Ma) in the above equation we use the 
Gaussian smoothing window with the radius determined 
by equation 13. The original PS formalism was based 
on the sharp k-space filter, but it has been shown sub- 
sequently by many authors that Gaussian filter works at 
least as well. The additional virtue of a Gaussian window 
is that it guarantees a convergent integral for <tq for any 
plausible power spectrum. Now, the only parameter left 
undetermined is 5 C , which will be fixed by comparing to 
N-body results. We find that the best overall fit of PS 
results to N-body results is obtained, if 

, _ fl.23- 0.05(1.0 -il z ) for A z = , 
c \1.23- 0.01(1.0 -Q z ) for fl z + A z = 1. {6 ) 

The results are shown in Figure 5. We see that PS 
fits N-body results quite well except for the Pk = k~ 2 
models (models 10,11,12,13). The PS results for all the 
= 1 models except the Pk = k~ 2 models appear to 
be somewhat above the N-body results at the low mass 
end (~ 5 x 10 14 /i _1 M©) and somewhat below the N-body 
results at the high mass end (~ 2.5 x 10 15 /i _1 M Q ). On 
the other hand, the PS results for all the Pk = k~ 2 mod- 
els are significantly above the N-body results. So, there 
is no room for further adjustments of 5 C to achieve better 
overall fits, at least for Gaussian smoothing windows. 

Note that S c (~ 1.23) is smaller than 1.67, which is 
in the expected direction because a smoother, Gaussian 
smoothing window is used here. 1.23 is also somewhat 
smaller than that given by Klypin et al. (1995), who give 
6 C = 1.40 for a Gaussian smoothing window in the context 
of damped Lyman alpha systems. But Klypin et al. also 
argue that S c could be as low as 1.3, were waves longer and 
shorter than those present in the simulation box included. 
We suspect that 5 C also depends on the shape of the power 
spectrum in a way that is analogous to the difference be- 
tween different smoothing windows: a steeper power spec- 
trum (i.e., n being smaller with Pk = k n ), which conspires 
to form a sharp k-space filter like that used in the original 
derivation of PS, requires a larger 8 C , while a flatter power 
spectrum requires a smaller S c . With this conjecture, the 
trend that applications of PS to smaller cosmic objects 
tend to require larger 6 C would have been predicted, since 
CDM-like spectra have a slowly bending shape which is 
flatter at small k (i.e., for larger systems) and steeper at 
large k (i.e., for smaller systems). This conjecture seems 
to be borne out in the subsequent analyses, as best sum- 
marized as in equation 36, where the dependence of as on 
T is consistent with the above hypothesis. This issue will 
be addressed elsewhere in more detail. 

While PS works well for CDM-like models for comput- 
ing halo mass functions, consistent with earlier works (Efs- 
tathiou & Rees 1988; WEF; ECF), it seems that GPM fits 
somewhat better the N-body results for CDM-like mod- 
els and also works for other models tested. An additional 
advantage is that GPM allows for a determination of the 
correlation function as well. Therefore, in subsequent cal- 
culations we will use GPM, if deemed appropriately appli- 
cable. 

3. VARIOUS FACTORS THAT AFFECT CLUSTER MASS 
FUNCTION 

It is worthwhile to understand what factors are rele- 
vant for the cluster mass function. We begin by show- 



ing the cluster mass functions for six variants of the stan- 
dard CDM model (Table 2 below), as indicated in panels 
(a,b,e,f,g,h) of Figure 6, at different normalization ampli- 
tudes (as)- We will return to Table 2 in §4 to discuss the 
various models in detail. The primordial power spectrum 
index is assumed to be n = 1 for the shown models in 
panels (a,b,e,f,g,h). Also shown as symbols are the obser- 
vations adopted from BC, and as three dashed curves are 
the fits to the symbols. The middle dashed curve is 

n fit (> M A JM x 10 _5 (M A /2.1 x 10 14 )" 1 

exp(-M A /2.1 x 10 14 ), (35) 

where Ufa is in /i 3 Mpc~ 3 and Ma is the cluster mass 
within the Abell radius in h~ x M©. This curve seems to 
represent the mean value of the observed mass function 
well (note that equation 35 is slightly different from the 
fitting formula in BC). The top and bottom dashed curves 
are 4n/ it (> Ma) and 0.25n/ it (> Ma)- It is difficult to 
estimate the errorbars of the observed mass function. The 
top and bottom dashed curves are intended to serve as 
2cr upper and lower bounds (in the vertical axis) of the 
observed mass function within the indicated mass range, 
which we deem to be conservative. Subsequent presenta- 
tions and explanations will follow this assertion. 

In all cases, we see that the cluster mass function be- 
comes progressively steeper at the high mass end as the 
amplitude of the density perturbations decreases. The 
physics behind this is simple to understand. As the am- 
plitude of fluctuations decreases, the required height of 
the density peaks for clusters with a given mass increases. 
Since the abundance of the high peaks at the very high 
end drops exponentially, the mass function steepens as 
the amplitude of fluctuations decreases. Note that there is 
only a narrow range in as where the model mass func- 
tion lies within the 2a limits in the mass range from 
4.8 x lO 14 ^" 1 M© to 1.2 x lO 15 ^ 1 M©. 

Although it is clear that the amplitude of the power 
spectrum sensitively determines the abundance of the clus- 
ters, as seen by comparing different curves within each 
panel, it is not yet clear how big the effect of f2o on the 
mass function is. In Figure 6, panel (c) is similar to panel 
(a) with only one change: f^o = 0.4 instead of £lo = 1.0. 
Note that the power spectrum used in panel (c) is identical 
to that used in panel (a). We see fio has a significant effect 
on the mass function. For example, the er 8 = 1.0 model 
in panel (c) has comparable mass function to the as = 0.6 
model in panel (a), but about two orders of magnitudes 
lower than the as — 1.0 model in panel (a). 

Next, we examine the effect of the shape of the power 
spectrum on the shape of the mass function. In Figure 6 
panel (d) is similar to panel (a) but with n = 0.75 instead 
of n — 1.0. Comparing panel (d) with panel (a) illustrates 
the sensitivity of the cluster mass function on the shape of 
power spectrum. The model shown in panel (d) has sub- 
stantially more power on large scales (100 — 300/i _1 Mpc) 
than the model shown in panel (a) . This difference results 
in flatter mass functions in panel (d) than in panel (a), es- 
pecially for the cases with lower as- A point that we would 
like to make here is that the cluster mass function depends 
not only on fi, A and as, but on Pk in a nontrivial, albeit 
relatively weak, way, especially at the high mass end of 
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the mass function. For example, at 1.2 x 10 15 ft. _1 Mq, the 
model shown in panel (d) has about a factor of four more 
clusters than the model shown in panel (a), with the only 
difference between the two models being a slight tilt of the 
power spectrum [n = 0.75 in (d) vsn = 1.0 in (a)]. 

Panels (e,g) and panels (f,h) show two pairs of low den- 
sity models, one being open and the other being flat with 
a cosmological constant but with a same Qq in each pair. 
One thing to note is that, other things being fixed, the 
cosmological constant does not make much difference in 
terms of the shape of the mass function. However, for a 
fixed <7g, the mass function in the model with a cosmo- 
logical constant is systematically lower than that of the 
model without a cosmological constant; the difference be- 
comes larger at lower amplitudes. 

Summarizing the above results we see that the factors 
that effect the cluster mass function in order of decreas- 
ing importance are ag, fioi Pk an d A . The ordering of 
the last two factors is somewhat more complex; the Pk 
factor is more significant at the high mass end, whereas 
the Ao factor affects rather uniformly across-the-board. It 
is therefore clear that the tightest constraint may be ob- 
tained for cr 8 for a given model with f2 , A and Pk being 
specified. 

4. NORMALIZING ALL CDM MODELS 

In §2.1 and 2.2 we have shown that the Gaussian peaks 
of cluster size indeed form clusters of galaxies and appro- 
priately identified peaks reproduce both cluster mass func- 
tion and cluster two-point correlation function accurately. 
We now use the observed zero redshift rich cluster abun- 
dance and the COBE observation to constrain six repre- 



sentative variants of the standard CDM model, which are 
of current interest. The models are listed in Table 2. 

The baryon densities for Models C,E are computed us- 
ing n b h 2 = 0.0125 (Walker et al. 1991), and for Models 
A,B,D,F using il b h? = 0.0193 (Buries & Tytler 1997). 
These choices of Q& for the models serve to maximize 
the viability of each model with respect to the observed 
gas fraction in X-ray clusters of galaxies (White et al. 
1993; Lubin et al. 1996; Danos & Pen 1998, p gas /ptot = 
(0.053 ± 0.004)ft,~ 3 / 2 ). The power spectrum transfer func- 
tions for all the models are computed using the CMBFAST 
code developed by Seljak and Zaldarriaga. The choice of 
the Hubble constant is made for each model such that 
each model is consistent with current measurements of 
the Hubble constant. It appears to be a consensus that 
Ho(obs) = 65 ± lOkm/s/Mpc can account for the dis- 
tribution of the current data from various measurements 
(see, e.g., Trimble 1997), except for those from Sunyaev- 
Zel'dovich observations (for a discussion of a reconciliation 
of this difference, see Cen 1998). Another consideration is 
that the age constraint from latest globular cluster ob- 
servations/interpretations (c.f., Salaris, Degl'Innocenti, & 
Weiss 1997) is not violated. 

The SCDM model (# 1 in Table 2) is the standard CDM 
with critical density. The mixed hot and cold dark matter 
model (HCDM) has critical density but with 20% of the 
mass in light massive neutrinos (two species of equal mass 
neutrinos arc assumed). The next two models, OCDM25 
and OCDM40 (# 3 and 4), are open models with matter 
density of f^o = 0.25 and 0.40, respectively. The last two 
models, ACDM25 and ACDM40 (# 5 and 6), are spatially 
flat models with the addition of a cosmological constant 



20 




(6e,f,g,h) 

Fig. 6. — Cluster mass functions for six models at several different normalization amplitudes of the power spectra. Also shown as symbols 
are the observations adopted from BC, and as three dashed curves are the fits to the symbols. The solid dots represent the cluster mass 
function from Abell cluster catalog, the open squares from the Edinburgh-Durham Cluster Catalog (Lumsden et al. 1992), and the open 
circles from the observed temperature function of Henry & Arnaud (1991). The middle dashed curve is computed by equation 35, nfa(> M), 
which represents the mean value of the observed mass function well. The top and bottom dashed curves are 4nfn(> M) and 0.25n^j t (> M), 
respectively. 



Table 2 
Six variants of CDM models 



Model Family 






A 


H 


fib 


tage (Gyrs) 


A 


0.936 


0.000 


0.000 


55 


0.064 


11.8 


B 


0.736 


0.200 


0.000 


55 


0.064 


11.8 


C 


0.220 


0.000 


0.000 


65 


0.030 


12.4 


D 


0.346 


0.000 


0.000 


60 


0.054 


12.7 


E 


0.220 


0.000 


0.750 


65 


0.030 


15.2 


F 


0.346 


0.000 


0.600 


60 


0.054 


14.5 
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but with otherwise similar parameters to the two open 
models. Note that for all models, a s and n are yet to 
be specified, as will be shown below, by normalizing each 
model to both COBE and zero redshift cluster abundance. 

GPM makes possible economically sample large (four di- 
mensional) parameter space spanned by the uncertainties 
in £lo, Ao, dg and Pk (the Hubble constant constitutes yet 
another dimension of uncertainty but its effect can be ab- 
sorbed into Pk in this case since the final quantities can all 
be expressed in units of h). The dependence on is sec- 
ondary and ignored here. If the cluster mass function of a 
model at any mass does not lie within the range delimited 
by the upper and lower dashed curves in Figure 6 in the 
mass range from 4.8 x lO 14 /^ 1 M to 1.2 x 10 15 /i _1 M Q , 
we conclude that the model is ruled out at 2a confidence 
level (this adoption of confidence level on the mass func- 
tion is somewhat crude but we think it is perhaps conser- 
vative). Thus, we can constrain the (Qo, A , as, Pk) four 
dimensional parameter space by the cluster mass function. 
The other, tight constraint currently available is from the 
COBE observations of the cosmic microwave background 
fluctuations on large scales. We will combine these two 
observations to constrain the models. 

With fio and Ao (and the composition of dark mat- 
ter and baryonic matter) being fixed, the power spectrum 
transfer function of a model becomes completely determin- 
istic. What is left to be specified is the primordial power 
spectrum index, n (assuming it is of a power law form). 
Figures (7a,b,c,d,e) show the constrained two-dimensional 
parameter space (as, n) for the six models listed in Table 2. 
The dotted hatched regions are the permitted space in this 
two-dimensional parameter plane whose cluster mass func- 
tions agree with observations at the 2a confidence level. 
The solid hatched regions are the constraints provided by 
COBE on large scales, where the spread in the vertical 
axis (as, ±14%, 2a) is due to the statistical uncertainty 
of the COBE measurements (Bunn & White 1997). We 
see that a quite tight constraint is obtained in the (n, as) 
plane. For the four spatially flat models (1,2,5,6) in Table 
2, two cases, with and without the tensor mode (gravi- 
tational wave) contribution to the CMB fluctuations on 
COBE scales, are considered. We assume the tensor to 
scalar ratio T/S = 7(1 — n) for n < 1 (Liddlc & Lyth 
1992; Davis et al. 1992; Crittenden et al. 1993; Stewart 
& Lyth 1993) and T/S = for n > 1 (Steinhardt 1997, 
private communications). For the remaining two models 
(open models; models 3,4 in Table 2) we only consider the 
case without the tensor mode contribution; it turns out 
to make no difference in the end whether tensor mode is 
considered or not, as the values of n in the allowed range 
are greater than unity. The permitted ranges of n, as and 
EP for all the models are tabulated in Table 3 to clearly 
show their distinctly different ranges for all the models. 
However, current observations by COBE on n (Gorski et 
al. 1996; Bennett et al. 1996; Hinshaw et al. 1996; Bunn 
& White 1997) and by large scale galaxy distribution on 
EP (Wright et al. 1992; or equivalently by T parameter, 
Maddox et al. 1990; Feldman, Kaiser, & Peacock 1994; 
Peacock, & Dodd 1994), as indicated as the last two rows 
in Table 3, while significantly constraining models at la 
confidence level, do not yet place 2a confidence level con- 
straint on this set of models. 

The values of n and as for the central model in each 



variant of the CDM model are shown as solid dots in six 
panels of Figure 7. To facilitate further examination of and 
intercomparison among the models, the parameters of the 
six central models are given in Table 4. These six models 
likely bracket all "viable" CDM models of current interest. 
This is an attempt to set the context for future discussions 
of CDM models that will share a common standard. 

5. cr 8 - RELATION 

Let us now steer off the main course for a moment and 
examine the as — relation that has been widely uti- 
lized (WEF; Viana & Liddlc 1996; ECF). We first note 
that it is very clear from Figure 7 and Table 3 that the 
shape of the power spectrum plays a significant role. While 
we find that the power spectrum shape dependence of a s 
for fl = 1 models (Models A,B in Table 3) is modest 
Aas ~ 0.03 — 0.1, in agreement with ECF, the power 
spectrum shape dependence of as for low density mod- 
els is substantial, being Aas — 0.17 — 0.30. Therefore, a 
more accurate relation of as — cannot be obtained until 
the shape of the power spectrum is more accurately fixed. 
The latter is currently unavailable, as indicated by the last 
two rows of Table 3. A tentative solution to circumvent 
this situation is to choose a shape of the power spectrum 
that is deemed to best fit observations of large-scale galaxy 
distribution. This approach was used by WEF, Viana & 
Liddle (1996) and ECF, which will also be adopted here. 

Figure 8 shows the constraint on the (Oo, as) plane, 
for two cases: T = 0.25 and A = (8a) and T = 0.25 
and fi + A = 1 (8b). The dotted region in each fig- 
ure represents the permitted regions of parameters where 
models have the cluster mass functions at z — consis- 
tent with the observations at 2a confidence level. The 
hatched region is the fit by ECF at 2a confidence level 
and the dotted curve is the fit by Viana & Liddle (1996). 
The dashed curve in (8b) is the fit by WEF. The solid 
curve is the best fit to the dotted region. The general 
agreements between various studies are very good. The 
small differences between the results of various studies 
perhaps reflect the differences in the procedures of fitting 
to observational data points, the small differences in the 
adopted observations being fit, and different theoretical 
approximations used for the models. Our best fit to the 
results of GPM to the observed mass function of BC, for 
r = 0.25, is (2a): a 8 (n Q ) = (0.52^;° 4 )fi °' 39 for A = 
and a 8 = (0.52t°;° 4 )«o - 43 for Q Q + A = 1. 

Given the sizable allowed range in T currently con- 
strained by observation (Peacock & Dodds 1994), T = 
0.23^0 03 (2a), it is useful to quantify the dependence 
of the above relationships on the shape of the power 
spectrum. Figure 8c shows the results and fits to them 
for T = 0.20 and T = 0.27, with the best fits be- 
ing (2a) <t 8 (Q ) = (0-485t^ 4 )^o"°' 39 for A = and 
a s = (0.485to ^)fi 43 for ft + A = 1 for V = 0.20, 
and cr 8 (Sl ) = (0.52±°;^)Slo °' 39 for A = and a s = 
(0.52±^°|)Oq °' 43 for Q + A = 1 for T = 0.27. The 
best fit of all the results for CDM-like models, allowing T 
to vary within the indicated range (0.20 — 0.27), is (2a): 

CO r w/[0-50 + 0.5(r- 0.23) ±0.05]«o 39 (A = 0) 

a8[ °' j ~li0.50 + 0.5(r-0.23)±0.05]fV - 43 (^o + A = 

(36) 
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Figure 7e 
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Fig. 7. — The constrained two-dimensional parameter space (a&,n) for the six models as tabulated in Table 2. The dotted hatched regions 
are the permitted space in this two-dimensional parameter plane whose cluster mass functions agree with observations at the 2<r confidence 
level. The solid hatched regions are the constraints provided by COBE on large scales, where the spread in the vertical axis (ag, ±14%, 2a) is 
due to the statistical uncertainty of the COBE measurements (Bunn & White 1997). For the four spatially flat models (models 1,2,5,6) two 
cases, with and without the tensor mode (gravitational wave) contribution to the CMB fluctuations on COBE scales are considered (see text). 
For the remaining two models (open models; models 3,4 of Table 2) we only consider the case without the gravitational wave contribution. 





Table 3 

Constraints on n, <t 8 and EP for the six variants of Table 2 



Model Family Tensor mode n range ag range EP range 



A 


yes 


0.72-0.82 


0.52- 


0.57 


1.01-1.05 


A 


no 


0.45 - 0.68 


0.49- 


0.52 


1.07- 1.20 


B 


yes 


0.79-0.92 


0.45- 


0.54 


1.44- 1.53 


B 


no 


0.61-0.84 


0.43- 


0.53 


1.49- 1.65 


C 


no 


1.39-1.49 


0.91 - 


1.14 


1.11 - 1.16 


D 


no 


0.98-1.23 


0.64- 


0.84 


1.25- 1.40 


E 


yes 


0.98-1.23 


0.81 - 


1.10 


1.25- 1.40 


E 


no 


0.96-1.23 


0.80- 


1.10 


1.25- 1.42 


F 


yes 


0.87-0.99 


0.65- 


0.82 


1.40- 1.47 


F 


no 


0.74-0.98 


0.62- 


0.81 


1.40- 1.56 


COBE(obs) 




1.20 ±0.60 (2a) 








LSS(obs) 










1.30 ±0.30 (2a) 



Table 4 

Six COBE and Cluster-Normalized CDM Models 



Model 


Ho 


n 


n c 




A 




^8 


tCDM 


55 


0.77 


0.936 


0.00 


0.0 


0.064 


0.55 


HCDM 


55 


0.88 


0.736 


0.20 


0.0 


0.064 


0.52 


OCDM25 


65 


1.47 


0.220 


0.00 


0.0 


0.030 


1.00 


OCDM40 


60 


1.15 


0.346 


0.00 


0.0 


0.054 


0.80 


ACDM25 


65 


1.10 


0.220 


0.00 


0.75 


0.030 


0.95 


ACDM40 


60 


0.96 


0.346 


0.00 


0.60 


0.054 


0.80 



Figure 8a 
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Figure 8b 





Fig. 8. — The constraint on the (f2o, its) plane, for two cases: V = 0.25 and Ao = (8a) and T = 0.25 and S7o + Ao = 1 (8b). The dotted 
region in each figure represents the permitted regions of parameters whose models have the cluster mass functions at z = consistent with 
what is observed at 2<r confidence level. The hatched region is the fit by Eke et al. (1996) at 2a confidence level and the short-dashed curve is 
the fit by Viana & Liddle (1996). The long-dashed curve is the fit by White et al. (1993). The solid curve is the best fit to the dotted region 
with the fitting formula being indicated in the figure. Figure 8c shows the results and fits to them for Y = 0.20 and V = 0.27. The best fit of 
the results, allowing V to vary within the indicated range, is presented in equation 36. 
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This us — 51o relation best summarizes a constraint on 
CDM-like models provided by the observations of zero red- 
shift rich cluster abundance. 

6. CONCLUSIONS 

We have developed and tested a method to compute the 
mass function and correlation function of peaks, based on 
the formalism for Gaussian density field. We should call 
this method Gaussian Peak Method (GPM). The essential 
new ingredient in this relatively old method is a simultane- 
ous determination of the smoothing window size (to select 
appropriate peaks) and the critical peak height of collapse. 
A large set of thirty-two N-body simulations are used to 
test the accuracy of the method and it is shown that the 
method is accurate for all the models tested, which cover 
the parameters space of interest spanned by Pj~, r2o, Ao 
and as- 

The GPM permits economical search of parameter 
space. We find that as — f^o relation is somewhat depen- 
dent upon the shape of the power spectrum. Normalizing 
CDM models to the observed local rich cluster abundance 
alone, allowing for the observed uncertainty in the shape 
parameter P(= Sl h)= (0.20 - 0.27), gives (2a): 

(Q r w/[0-50 + 0.5(r- 0.23) ±0.05]^ 39 (A = 0) 
as( °' ' 110.50 + 0.5(r- 0.23) iO.05]^ 43 (O + A = 

(37) 



Matching both COBE on very large scales and the abun- 
dance of local rich clusters of galaxies fixes both the shape 
(n) and amplitude of the power spectrum (eg) of any 
model to about 10% accuracy. Consequently, all mod- 
els become almost completely deterministic. A set of six 
CDM models (including cold plus hot dark matter model) 
(Table 4), likely bracketing all potentially interesting mod- 
els, is advertized together. This should be viewed as an 
attempt to set the context for future discussions on a set 
of standardized models to facilitate comparison of results 
between workers in the field. 
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